#=
测试哈密顿量构造
=#


include("../src/edmary.jl")

using ..edmary

using SparseArrays
using LinearAlgebra
using Arpack: eigs
using Dates
using JLD2


function saveCSC(fh, name)
    if isfile(name)
        rm(name)
    end
    handle = jldopen(name, "w")
    write(handle, "CSC", fh)
    close(handle)
end


function loadCSC(name)
    handle = jldopen(name, "r")
    csc = read(handle, "CSC")
    close(handle)
    return csc
end


function run()
    L = 4
    npart = 3
    nsite = L^2
    U = 3.0
    t1 = -1.0
    t2 = -1.5
    name = "L$(L)npart$(npart)nsite$(nsite)U$(U)t1$(t1)t2$(t2).jld"
    sta, invsta = get_state_array(npart, nsite)
    println(length(sta))
    @show Dates.now()
    bonds = bonds_ckb(L, t1, t2)
    if isfile(name)
        fh = loadCSC(name)
    else
        fh = full_hmltCSC(npart, nsite, sta, invsta, bonds, U)
        saveCSC(fh, name)
    end
    @show Dates.now()
    #println(Matrix(fh))
    #
    #nurand = rand(length(sta))
    #nutmp = zeros(length(sta))
    #apply_hmlt!(nutmp, nurand, nsite, sta, invsta, bonds, U)
    #println(nutmp)
    #println(fh*nurand)
    #@show Dates.now()
    #@assert all(isapprox.(nutmp, fh*nurand, atol=1e-10))
    #@show Dates.now()
    #
    v0r = rand(length(sta))
    v0r = v0r / norm(v0r)
    eval, evec = eigs(fh; nev=2, which=:SR, tol=0.0, maxiter=100, sigma=nothing, ritzvec=true, v0=v0r)
    println(eval)
    @show Dates.now()
    #
    coef1 = sqrt(3)/2
    coef2 = complex(0.0, 1/2)
    ppbonds = ppbonds_ckb(L,
    Dict("+x"=>coef1, "+y"=>-coef1, "-x"=>coef1, "-y"=>-coef1, "+x+y"=>coef2, "-x-y"=>coef2),
    Dict("+x"=>coef1, "+y"=>-coef1, "-x"=>coef1, "-y"=>-coef1, "+x-y"=>coef2, "-x+y"=>coef2),
    )
    corr = ppbuild_from_order_param(ppbonds, [1, 1+L])
    gf = zeros(nsite, nsite)
    #
    gfu = zeros(nsite, nsite)
    gfd = zeros(nsite, nsite)
    for i in Base.OneTo(nsite); for j in Base.OneTo(nsite)
        gfciju = gfc_at_ij(i, j, :up, npart, nsite, sta, invsta)
        gndv = @view(evec[:, 1])
        pvec = gfciju*gndv
        gfu[i, j] = adjoint(gndv)*pvec
    end; end
    #
    for i in Base.OneTo(nsite); for j in Base.OneTo(nsite)
        gfcijd = gfc_at_ij(i, j, :dn, npart, nsite, sta, invsta)
        gndv = @view(evec[:, 1])
        pvec = gfcijd*gndv
        gfd[i, j] = adjoint(gndv)*pvec
    end; end
    @assert all(isapprox.(gfu, gfd, atol=1e-8))
    gf = gfu
    println(gf)
    unc = uncorr_from_ppbuild(corr, gf)
    println(unc)
end


run()
